data <- readRDS("data/raw/very_low_birthweight.RDS")
# Добавление id
data <- data %>%
mutate(id = row_number())
# Перевод категориальных переменных в факторы
data <- data %>% mutate(
across(c(id, race, inout, delivery, pvh, ivh, ipe, sex), ~ as.factor(.x)))
cleaned_data <- data %>% select(where(~ sum(is.na(.)) <= 100)) %>%
filter(across(everything(), ~ !is.na(.)))
## Warning: Using `across()` in `filter()` was deprecated in dplyr 1.0.8.
## ℹ Please use `if_any()` or `if_all()` instead.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
outliers <- function(x) {
(x < quantile(x, 0.25) - 3*IQR(x)) |
(x > quantile(x, 0.75) + 3*IQR(x))
}
cleaned_data %>% filter(across(where(is.numeric), ~ !outliers(.))) %>% glimpse()
## Warning: Using `across()` in `filter()` was deprecated in dplyr 1.0.8.
## ℹ Please use `if_any()` or `if_all()` instead.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Rows: 277
## Columns: 20
## $ birth <dbl> 81.514, 81.558, 81.610, 81.624, 81.626, 81.689, 81.697, 81.70…
## $ exit <dbl> 81.539, 81.667, 81.697, 81.700, 81.730, 81.878, 81.952, 81.82…
## $ hospstay <int> 9, 40, 32, 28, 38, 69, 93, 44, 44, 70, 85, 58, 75, 35, 34, 46…
## $ lowph <dbl> 7.250000, 7.250000, 7.320000, 7.160000, 7.039997, 7.419998, 7…
## $ pltct <int> 244, 182, 282, 153, 229, 361, 255, 186, 134, 229, 68, 174, 17…
## $ race <fct> white, black, black, black, white, white, black, white, white…
## $ bwt <int> 1370, 1480, 1255, 1350, 1310, 1180, 770, 1490, 1000, 1120, 74…
## $ gest <dbl> 32.0, 32.0, 29.5, 34.0, 32.0, 28.0, 26.0, 33.0, 28.0, 29.0, 2…
## $ inout <fct> born at Duke, born at Duke, born at Duke, born at Duke, born …
## $ twn <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ delivery <fct> abdominal, vaginal, vaginal, abdominal, vaginal, abdominal, v…
## $ apg1 <int> 7, 8, 9, 4, 6, 6, 4, 8, 5, 9, 9, 9, 1, 8, 4, 6, 9, 6, 8, 9, 8…
## $ vent <int> 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1…
## $ pneumo <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ pda <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ cld <int> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0…
## $ year <dbl> 81.51471, 81.55847, 81.61053, 81.62421, 81.62695, 81.68988, 8…
## $ sex <fct> female, male, female, female, male, male, male, male, female,…
## $ dead <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ id <fct> 2, 4, 7, 10, 11, 14, 16, 17, 21, 22, 23, 25, 31, 35, 36, 41, …
data_long <- cleaned_data %>%
select(where(is.numeric)) %>%
pivot_longer(everything(), names_to = "variable", values_to = "value")
ggplot(data_long, aes(x = value)) +
geom_density(fill = "blue") +
facet_wrap(~ variable, scales = "free")+
theme_bw()
# Графики плотности для переменных bwt и gest с окраской по inout
data_long <- cleaned_data %>%
pivot_longer(cols = c(bwt, gest), names_to = "variable", values_to = "value")
ggplot(data_long, aes(x = value, fill = inout)) +
geom_density(alpha = 0.5) +
facet_wrap(~ variable, scales = "free", nrow=2)+
theme_bw()
# Cравнение значений колонки ‘lowph’ между группами в переменной
inout
H0: среднее lowph в группе born at Duke равно среднему в группе transported
H1: средние не равны
Большие выборки, распределение приближается к нормальному по ЦПТ, используем двусторонний t тест c поправкой Уэлча для неизвестных дисперсий (уровень значимости установим 0.05)
t = 5.5731, df = 111.03, p-value = 1.77e-07
p-value<0.05
95 percent confidence interval: 0.05762847 0.1212193
0 не входит в ДИ
Следовательно, отклоняем нулевую гипотезу в пользу альтернативной
Как бы вы интерпретировали результат, если бы знали, что более низкое значение lowph ассоциировано с более низкой выживаемостью?
В таком случае группа transported имеет более низкую выживаемость по сравнению с группой born at Duke
test <- cleaned_data %>% t_test(lowph ~ inout, detailed=TRUE) %>%
add_significance() %>%
mutate(y.position = max(cleaned_data$lowph, na.rm = TRUE) + 0.5)
## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "born at Duke" - "transported", or divided in the order "born at Duke" /
## "transported" for ratio-based statistics. To specify this order yourself,
## supply `order = c("born at Duke", "transported")`.
print(test)
## # A tibble: 1 × 8
## statistic t_df p_value alternative estimate lower_ci upper_ci y.position
## <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 5.57 111. 0.000000177 two.sided 0.0894 0.0576 0.121 8.05
ggplot(cleaned_data, aes(x=inout, y=lowph, fill=inout))+
geom_boxplot()+
#stat_pvalue_manual(test, label="p-value")+ не работает с этой строкой, в чем причина?
theme_bw()
num_data <- cleaned_data %>% select(where(is.numeric), -c(birth,year, exit)) %>% glimpse()
## Rows: 531
## Columns: 12
## $ hospstay <int> 9, 40, 2, 32, 28, 38, 62, 69, 1, 93, 44, 66, 65, 44, 70, 85, …
## $ lowph <dbl> 7.250000, 7.250000, 6.969997, 7.320000, 7.160000, 7.039997, 7…
## $ pltct <int> 244, 182, 54, 282, 153, 229, 182, 361, 378, 255, 186, 260, 18…
## $ bwt <int> 1370, 1480, 925, 1255, 1350, 1310, 1110, 1180, 970, 770, 1490…
## $ gest <dbl> 32.0, 32.0, 28.0, 29.5, 34.0, 32.0, 28.0, 28.0, 28.0, 26.0, 3…
## $ twn <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1…
## $ apg1 <int> 7, 8, 5, 9, 4, 6, 6, 6, 2, 4, 8, 1, 8, 5, 9, 9, 9, 6, 2, 1, 6…
## $ vent <int> 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0…
## $ pneumo <int> 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0…
## $ pda <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ cld <int> 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0…
## $ dead <int> 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0…
data_cor <- cor(num_data)
data_cor %>% corr.test(method = "spearman") %>%
print(short=FALSE)
## Call:corr.test(x = ., method = "spearman")
## Correlation matrix
## hospstay lowph pltct bwt gest twn apg1 vent pneumo pda cld
## hospstay 1.00 -0.76 -0.71 -0.87 -0.83 -0.61 -0.76 0.84 0.64 0.83 0.95
## lowph -0.76 1.00 0.83 0.85 0.90 0.60 0.84 -0.97 -0.94 -0.85 -0.85
## pltct -0.71 0.83 1.00 0.78 0.75 0.41 0.91 -0.85 -0.83 -0.82 -0.80
## bwt -0.87 0.85 0.78 1.00 0.97 0.80 0.92 -0.90 -0.74 -0.92 -0.95
## gest -0.83 0.90 0.75 0.97 1.00 0.83 0.87 -0.94 -0.80 -0.94 -0.94
## twn -0.61 0.60 0.41 0.80 0.83 1.00 0.65 -0.64 -0.39 -0.65 -0.71
## apg1 -0.76 0.84 0.91 0.92 0.87 0.65 1.00 -0.86 -0.76 -0.83 -0.86
## vent 0.84 -0.97 -0.85 -0.90 -0.94 -0.64 -0.86 1.00 0.88 0.91 0.93
## pneumo 0.64 -0.94 -0.83 -0.74 -0.80 -0.39 -0.76 0.88 1.00 0.83 0.73
## pda 0.83 -0.85 -0.82 -0.92 -0.94 -0.65 -0.83 0.91 0.83 1.00 0.94
## cld 0.95 -0.85 -0.80 -0.95 -0.94 -0.71 -0.86 0.93 0.73 0.94 1.00
## dead 0.69 -0.97 -0.79 -0.82 -0.87 -0.55 -0.80 0.90 0.97 0.85 0.78
## dead
## hospstay 0.69
## lowph -0.97
## pltct -0.79
## bwt -0.82
## gest -0.87
## twn -0.55
## apg1 -0.80
## vent 0.90
## pneumo 0.97
## pda 0.85
## cld 0.78
## dead 1.00
## Sample Size
## [1] 12
## Probability values (Entries above the diagonal are adjusted for multiple tests.)
## hospstay lowph pltct bwt gest twn apg1 vent pneumo pda cld dead
## hospstay 0.00 0.07 0.11 0.01 0.03 0.20 0.08 0.02 0.20 0.03 0.00 0.13
## lowph 0.00 0.00 0.03 0.02 0.00 0.20 0.02 0.00 0.00 0.02 0.02 0.00
## pltct 0.01 0.00 0.00 0.06 0.08 0.37 0.00 0.02 0.03 0.03 0.05 0.05
## bwt 0.00 0.00 0.00 0.00 0.00 0.05 0.00 0.00 0.08 0.00 0.00 0.03
## gest 0.00 0.00 0.01 0.00 0.00 0.03 0.01 0.00 0.05 0.00 0.00 0.01
## twn 0.04 0.04 0.18 0.00 0.00 0.00 0.20 0.20 0.37 0.20 0.11 0.20
## apg1 0.00 0.00 0.00 0.00 0.00 0.02 0.00 0.01 0.08 0.03 0.01 0.05
## vent 0.00 0.00 0.00 0.00 0.00 0.02 0.00 0.00 0.01 0.00 0.00 0.00
## pneumo 0.03 0.00 0.00 0.01 0.00 0.21 0.00 0.00 0.00 0.03 0.10 0.00
## pda 0.00 0.00 0.00 0.00 0.00 0.02 0.00 0.00 0.00 0.00 0.00 0.02
## cld 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.01 0.00 0.00 0.06
## dead 0.01 0.00 0.00 0.00 0.00 0.06 0.00 0.00 0.00 0.00 0.00 0.00
##
## Confidence intervals based upon normal theory. To get bootstrapped values, try cor.ci
## raw.lower raw.r raw.upper raw.p lower.adj upper.adj
## hspst-lowph -0.93 -0.76 -0.33 0.00 -0.96 0.00
## hspst-pltct -0.91 -0.71 -0.22 0.01 -0.95 0.07
## hspst-bwt -0.96 -0.87 -0.58 0.00 -0.98 -0.23
## hspst-gest -0.95 -0.83 -0.48 0.00 -0.98 -0.13
## hspst-twn -0.88 -0.61 -0.05 0.04 -0.92 0.15
## hspst-apg1 -0.93 -0.76 -0.32 0.00 -0.96 0.00
## hspst-vent 0.51 0.84 0.95 0.00 0.15 0.98
## hspst-pneum 0.10 0.64 0.89 0.03 -0.13 0.93
## hspst-pda 0.49 0.83 0.95 0.00 0.14 0.98
## hspst-cld 0.83 0.95 0.99 0.00 0.62 0.99
## hspst-dead 0.20 0.69 0.91 0.01 -0.08 0.95
## lowph-pltct 0.49 0.83 0.95 0.00 0.14 0.98
## lowph-bwt 0.55 0.85 0.96 0.00 0.19 0.98
## lowph-gest 0.68 0.90 0.97 0.00 0.37 0.99
## lowph-twn 0.04 0.60 0.87 0.04 -0.14 0.91
## lowph-apg1 0.51 0.84 0.95 0.00 0.15 0.98
## lowph-vent -0.99 -0.97 -0.88 0.00 -1.00 -0.71
## lowph-pneum -0.98 -0.94 -0.81 0.00 -0.99 -0.58
## lowph-pda -0.96 -0.85 -0.55 0.00 -0.98 -0.19
## lowph-cld -0.96 -0.85 -0.53 0.00 -0.98 -0.17
## lowph-dead -0.99 -0.97 -0.88 0.00 -1.00 -0.71
## pltct-bwt 0.36 0.78 0.93 0.00 0.03 0.97
## pltct-gest 0.31 0.75 0.92 0.01 -0.01 0.96
## pltct-twn -0.21 0.41 0.80 0.18 -0.30 0.83
## pltct-apg1 0.70 0.91 0.97 0.00 0.40 0.99
## pltct-vent -0.96 -0.85 -0.55 0.00 -0.98 -0.19
## pltct-pneum -0.95 -0.83 -0.48 0.00 -0.98 -0.13
## pltct-pda -0.95 -0.82 -0.46 0.00 -0.98 -0.12
## pltct-cld -0.94 -0.80 -0.41 0.00 -0.97 -0.07
## pltct-dead -0.94 -0.79 -0.40 0.00 -0.97 -0.06
## bwt-gest 0.90 0.97 0.99 0.00 0.76 1.00
## bwt-twn 0.41 0.80 0.94 0.00 0.07 0.97
## bwt-apg1 0.74 0.92 0.98 0.00 0.47 0.99
## bwt-vent -0.97 -0.90 -0.66 0.00 -0.99 -0.34
## bwt-pneum -0.92 -0.74 -0.29 0.01 -0.96 0.02
## bwt-pda -0.98 -0.92 -0.74 0.00 -0.99 -0.47
## bwt-cld -0.99 -0.95 -0.83 0.00 -0.99 -0.62
## bwt-dead -0.95 -0.82 -0.46 0.00 -0.98 -0.11
## gest-twn 0.48 0.83 0.95 0.00 0.12 0.98
## gest-apg1 0.58 0.87 0.96 0.00 0.23 0.98
## gest-vent -0.98 -0.94 -0.79 0.00 -0.99 -0.54
## gest-pneum -0.94 -0.80 -0.41 0.00 -0.97 -0.06
## gest-pda -0.98 -0.94 -0.79 0.00 -0.99 -0.54
## gest-cld -0.98 -0.94 -0.79 0.00 -0.99 -0.54
## gest-dead -0.96 -0.87 -0.58 0.00 -0.98 -0.23
## twn-apg1 0.12 0.65 0.89 0.02 -0.13 0.93
## twn-vent -0.89 -0.64 -0.11 0.02 -0.93 0.13
## twn-pneum -0.79 -0.39 0.24 0.21 -0.79 0.24
## twn-pda -0.89 -0.65 -0.12 0.02 -0.94 0.15
## twn-cld -0.91 -0.71 -0.24 0.01 -0.95 0.06
## twn-dead -0.86 -0.55 0.03 0.06 -0.89 0.17
## apg1-vent -0.96 -0.86 -0.57 0.00 -0.98 -0.21
## apg1-pneum -0.93 -0.76 -0.32 0.00 -0.96 0.01
## apg1-pda -0.95 -0.83 -0.49 0.00 -0.98 -0.13
## apg1-cld -0.96 -0.86 -0.57 0.00 -0.98 -0.21
## apg1-dead -0.94 -0.80 -0.41 0.00 -0.97 -0.06
## vent-pneum 0.62 0.88 0.97 0.00 0.28 0.99
## vent-pda 0.70 0.91 0.97 0.00 0.40 0.99
## vent-cld 0.76 0.93 0.98 0.00 0.50 0.99
## vent-dead 0.66 0.90 0.97 0.00 0.34 0.99
## pneum-pda 0.48 0.83 0.95 0.00 0.12 0.98
## pneum-cld 0.26 0.73 0.92 0.01 -0.04 0.96
## pneum-dead 0.88 0.97 0.99 0.00 0.71 1.00
## pda-cld 0.79 0.94 0.98 0.00 0.54 0.99
## pda-dead 0.53 0.85 0.96 0.00 0.17 0.98
## cld-dead 0.36 0.78 0.93 0.00 0.03 0.97
corrplot(data_cor, method = 'number',insig = "blank")
data_cor %>%
network_plot(min_cor = .0)
# Иерархическая кластеризация
# Стандартизация
data_scaled <- scale(data_cor)
# Матрица дистанций
data_dist <- dist(data_scaled, method = "euclidean")
# Дендрограмма кластеров
data_dist.hc <- hclust(d = data_dist, method = "ward.D2")
# Визуализация
fviz_dend(data_dist.hc,
k=2,
cex = 0.6,
color_labels_by_k = TRUE,
rect = TRUE)
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
## Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Определение оптимального количества кластеров методом силуэт и потом оптимизация кода выше
fviz_nbclust(data_scaled, kmeans, method = "silhouette")
# Heatmap+кластеризация
Вывод: в данных имеется два кластера: 1(pneumo, dead, hospstay, pda, vent, cld) и 2 (lowph, bwt, gest, pltct, apg1, twn). Между кластерами наблюдается преимущественно отрицательная корреляция.
pheatmap(data_scaled,
show_rownames = TRUE,
clustering_distance_rows = data_dist,
clustering_method = "ward.D2",
cutree_rows = 5,
cutree_cols = length(colnames(data_scaled)),
angle_col = 45,
main = "Dendrograms for clustering rows and columns with heatmap")
# PCA
Вывод: первая компонента описывает 65% вариаций, а первые две компоненты описывают 75.45% вариаций данных (следовательно, в этих данных много скоррелированных между собой переменных). 1 группа (pneumo, dead, hospstay, pda, vent, cld) и 2 (lowph, bwt, gest, pltct, apg1, twn) имеют отрицательную корреляцию между группами, но при этом они скоррелированы между собой внутри групп.
# Используем стандартизованные значения для PCA, так как шкалы у переменных разные
data_pca <- prcomp(num_data, scale = T)
summary(data_pca)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.9667 1.1625 1.07954 0.99768 0.90725 0.88700 0.83410
## Proportion of Variance 0.3223 0.1126 0.09712 0.08295 0.06859 0.06556 0.05798
## Cumulative Proportion 0.3223 0.4350 0.53208 0.61503 0.68362 0.74918 0.80716
## PC8 PC9 PC10 PC11 PC12
## Standard deviation 0.81596 0.76417 0.66649 0.60125 0.50858
## Proportion of Variance 0.05548 0.04866 0.03702 0.03013 0.02155
## Cumulative Proportion 0.86264 0.91130 0.94832 0.97845 1.00000
fviz_pca_var(data_pca, col.var = "contrib")
# biplots
biplot <- ggbiplot(data_pca,
groups = as.factor(cleaned_data$dead),
scale=0, alpha = 0.5) +
theme_bw()
biplot
data_pca$id <-cleaned_data$id
ggplotly(biplot) %>%
style(text = paste("ID:", data_pca$id)) # как сделать без перекрытия переменных?
Колонку dead использовать некорректно, так как она не учитывает время жизни (exit-birth), что необходимо для модели выживаемости
Вывод: Два кластера (голубой и красный в углу) отличаются от остальных по каким-то признакам, при этом часть живых пациентов имеют характеристики, близкие к умершим. То же самое видно на PCA (умершие имеют позитивную корреляцию с одним кластером переменных и обратную корреляцию с другим кластером переменных, но в то же время часть живых пациентов тоже имеют такую зависимость).
umap_prep <- recipe(~., data = num_data) %>%
step_normalize(all_predictors()) %>%
step_umap(all_predictors()) %>%
prep() %>%
juice()
umap_prep %>%
ggplot(aes(UMAP1, UMAP2)) +
geom_point(aes(color = as.character(num_data$dead)),
alpha = 0.7, size = 1) +
labs(color = NULL)
# Изменим n_neighbors и min_dist
Вывод: появилось 2 голубых кластера, данные выглядят более однородными
umap_prep <- recipe(~., data = num_data) %>%
step_normalize(all_predictors()) %>%
step_umap(all_predictors(),neighbors=5, min_dist=0.5) %>% # по умолчанию было 15 и 0.01 соответственно
prep() %>%
juice()
umap_prep %>%
ggplot(aes(UMAP1, UMAP2)) +
geom_point(aes(color = as.character(num_data$dead)),
alpha = 0.7, size = 1) +
labs(color = NULL)